library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(raster)
## Loading required package: sp
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(stringr) # for working with strings (pattern matching)
library(tidyr) # for unite() and separate()
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:raster':
##
## extract
library(spData)
world
## Simple feature collection with 177 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS: WGS 84
## # A tibble: 177 x 11
## iso_a2 name_long continent region_un subregion type area_km2 pop lifeExp
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 FJ Fiji Oceania Oceania Melanesia Sove… 1.93e4 8.86e5 70.0
## 2 TZ Tanzania Africa Africa Eastern … Sove… 9.33e5 5.22e7 64.2
## 3 EH Western … Africa Africa Northern… Inde… 9.63e4 NA NA
## 4 CA Canada North Am… Americas Northern… Sove… 1.00e7 3.55e7 82.0
## 5 US United S… North Am… Americas Northern… Coun… 9.51e6 3.19e8 78.8
## 6 KZ Kazakhst… Asia Asia Central … Sove… 2.73e6 1.73e7 71.6
## 7 UZ Uzbekist… Asia Asia Central … Sove… 4.61e5 3.08e7 71.0
## 8 PG Papua Ne… Oceania Oceania Melanesia Sove… 4.65e5 7.76e6 65.2
## 9 ID Indonesia Asia Asia South-Ea… Sove… 1.82e6 2.55e8 68.9
## 10 AR Argentina South Am… Americas South Am… Sove… 2.78e6 4.30e7 76.3
## # … with 167 more rows, and 2 more variables: gdpPercap <dbl>,
## # geom <MULTIPOLYGON [°]>
coffee_data
## # A tibble: 47 x 3
## name_long coffee_production_2016 coffee_production_2017
## <chr> <int> <int>
## 1 Angola NA NA
## 2 Bolivia 3 4
## 3 Brazil 3277 2786
## 4 Burundi 37 38
## 5 Cameroon 8 6
## 6 Central African Republic NA NA
## 7 Congo, Dem. Rep. of 4 12
## 8 Colombia 1330 1169
## 9 Costa Rica 28 32
## 10 Côte d'Ivoire 114 130
## # … with 37 more rows
world_coffee = left_join(world, coffee_data, by = "name_long") #2つのdata.frameの共通点を元にくっつける。byで明示
#> Joining, by = "name_long"
class(world_coffee)
## [1] "sf" "tbl_df" "tbl" "data.frame"
#> [1] "sf" "tbl_df" "tbl" "data.frame"
world_coffee
## Simple feature collection with 177 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS: WGS 84
## # A tibble: 177 x 13
## iso_a2 name_long continent region_un subregion type area_km2 pop lifeExp
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 FJ Fiji Oceania Oceania Melanesia Sove… 1.93e4 8.86e5 70.0
## 2 TZ Tanzania Africa Africa Eastern … Sove… 9.33e5 5.22e7 64.2
## 3 EH Western … Africa Africa Northern… Inde… 9.63e4 NA NA
## 4 CA Canada North Am… Americas Northern… Sove… 1.00e7 3.55e7 82.0
## 5 US United S… North Am… Americas Northern… Coun… 9.51e6 3.19e8 78.8
## 6 KZ Kazakhst… Asia Asia Central … Sove… 2.73e6 1.73e7 71.6
## 7 UZ Uzbekist… Asia Asia Central … Sove… 4.61e5 3.08e7 71.0
## 8 PG Papua Ne… Oceania Oceania Melanesia Sove… 4.65e5 7.76e6 65.2
## 9 ID Indonesia Asia Asia South-Ea… Sove… 1.82e6 2.55e8 68.9
## 10 AR Argentina South Am… Americas South Am… Sove… 2.78e6 4.30e7 76.3
## # … with 167 more rows, and 4 more variables: gdpPercap <dbl>,
## # geom <MULTIPOLYGON [°]>, coffee_production_2016 <int>,
## # coffee_production_2017 <int>
plot(world_coffee["coffee_production_2017"])
world_coffee_inner <- inner_join(world, coffee_data) #inner_join -> left_joinでNAになった要素を弾く
## Joining, by = "name_long"
world_coffee_inner
## Simple feature collection with 45 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -117.1278 ymin: -33.76838 xmax: 156.02 ymax: 35.49401
## Geodetic CRS: WGS 84
## # A tibble: 45 x 13
## iso_a2 name_long continent region_un subregion type area_km2 pop lifeExp
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 TZ Tanzania Africa Africa Eastern … Sove… 932746. 5.22e7 64.2
## 2 PG Papua New… Oceania Oceania Melanesia Sove… 464520. 7.76e6 65.2
## 3 ID Indonesia Asia Asia South-Ea… Sove… 1819251. 2.55e8 68.9
## 4 KE Kenya Africa Africa Eastern … Sove… 590837. 4.60e7 66.2
## 5 DO Dominican… North Am… Americas Caribbean Sove… 48158. 1.04e7 73.5
## 6 TL Timor-Les… Asia Asia South-Ea… Sove… 14715. 1.21e6 68.3
## 7 MX Mexico North Am… Americas Central … Sove… 1969480. 1.24e8 76.8
## 8 BR Brazil South Am… Americas South Am… Sove… 8508557. 2.04e8 75.0
## 9 BO Bolivia South Am… Americas South Am… Sove… 1085270. 1.06e7 68.4
## 10 PE Peru South Am… Americas South Am… Sove… 1309700. 3.10e7 74.5
## # … with 35 more rows, and 4 more variables: gdpPercap <dbl>,
## # geom <MULTIPOLYGON [°]>, coffee_production_2016 <int>,
## # coffee_production_2017 <int>
plot(world_coffee_inner["coffee_production_2017"])
world_new = world # do not overwrite our original data
world_new$pop_dens = world_new$pop / world_new$area_km2 #計算されたpop_densをworld_newに追加
world_new
## Simple feature collection with 177 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS: WGS 84
## # A tibble: 177 x 12
## iso_a2 name_long continent region_un subregion type area_km2 pop lifeExp
## * <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 FJ Fiji Oceania Oceania Melanesia Sove… 1.93e4 8.86e5 70.0
## 2 TZ Tanzania Africa Africa Eastern … Sove… 9.33e5 5.22e7 64.2
## 3 EH Western … Africa Africa Northern… Inde… 9.63e4 NA NA
## 4 CA Canada North Am… Americas Northern… Sove… 1.00e7 3.55e7 82.0
## 5 US United S… North Am… Americas Northern… Coun… 9.51e6 3.19e8 78.8
## 6 KZ Kazakhst… Asia Asia Central … Sove… 2.73e6 1.73e7 71.6
## 7 UZ Uzbekist… Asia Asia Central … Sove… 4.61e5 3.08e7 71.0
## 8 PG Papua Ne… Oceania Oceania Melanesia Sove… 4.65e5 7.76e6 65.2
## 9 ID Indonesia Asia Asia South-Ea… Sove… 1.82e6 2.55e8 68.9
## 10 AR Argentina South Am… Americas South Am… Sove… 2.78e6 4.30e7 76.3
## # … with 167 more rows, and 3 more variables: gdpPercap <dbl>,
## # geom <MULTIPOLYGON [°]>, pop_dens <dbl>
world_new2 <- world %>%
mutate(pop_dens = pop/area_km2)
world_new2 #上と同じ
## Simple feature collection with 177 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS: WGS 84
## # A tibble: 177 x 12
## iso_a2 name_long continent region_un subregion type area_km2 pop lifeExp
## * <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 FJ Fiji Oceania Oceania Melanesia Sove… 1.93e4 8.86e5 70.0
## 2 TZ Tanzania Africa Africa Eastern … Sove… 9.33e5 5.22e7 64.2
## 3 EH Western … Africa Africa Northern… Inde… 9.63e4 NA NA
## 4 CA Canada North Am… Americas Northern… Sove… 1.00e7 3.55e7 82.0
## 5 US United S… North Am… Americas Northern… Coun… 9.51e6 3.19e8 78.8
## 6 KZ Kazakhst… Asia Asia Central … Sove… 2.73e6 1.73e7 71.6
## 7 UZ Uzbekist… Asia Asia Central … Sove… 4.61e5 3.08e7 71.0
## 8 PG Papua Ne… Oceania Oceania Melanesia Sove… 4.65e5 7.76e6 65.2
## 9 ID Indonesia Asia Asia South-Ea… Sove… 1.82e6 2.55e8 68.9
## 10 AR Argentina South Am… Americas South Am… Sove… 2.78e6 4.30e7 76.3
## # … with 167 more rows, and 3 more variables: gdpPercap <dbl>,
## # geom <MULTIPOLYGON [°]>, pop_dens <dbl>
world %>% transmute(pop_dens = pop/area_km2) #計算されたものだけ残る
## Simple feature collection with 177 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## Geodetic CRS: WGS 84
## # A tibble: 177 x 2
## pop_dens geom
## * <dbl> <MULTIPOLYGON [°]>
## 1 45.9 (((180 -16.06713, 180 -16.55522, 179.3641 -16.80135, 178.7251 -17.0…
## 2 56.0 (((33.90371 -0.95, 34.07262 -1.05982, 37.69869 -3.09699, 37.7669 -3…
## 3 NA (((-8.66559 27.65643, -8.665124 27.58948, -8.6844 27.39574, -8.6872…
## 4 3.54 (((-122.84 49, -122.9742 49.00254, -124.9102 49.98456, -125.6246 50…
## 5 33.5 (((-122.84 49, -120 49, -117.0312 49, -116.0482 49, -113 49, -110.0…
## 6 6.33 (((87.35997 49.21498, 86.59878 48.54918, 85.76823 48.45575, 85.7204…
## 7 66.7 (((55.96819 41.30864, 55.92892 44.99586, 58.50313 45.5868, 58.68999…
## 8 16.7 (((141.0002 -2.600151, 142.7352 -3.289153, 144.584 -3.861418, 145.2…
## 9 140. (((141.0002 -2.600151, 141.0171 -5.859022, 141.0339 -9.117893, 140.…
## 10 15.4 (((-68.63401 -52.63637, -68.25 -53.1, -67.75 -53.85, -66.45 -54.45,…
## # … with 167 more rows
world %>% transmute(pop_dens = pop/area_km2) %>% st_drop_geometry()
## # A tibble: 177 x 1
## pop_dens
## * <dbl>
## 1 45.9
## 2 56.0
## 3 NA
## 4 3.54
## 5 33.5
## 6 6.33
## 7 66.7
## 8 16.7
## 9 140.
## 10 15.4
## # … with 167 more rows
#raster画像を自分で作成
elev = raster(nrows = 6, ncols = 6, res = 0.5,
xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
vals = 1:36)
plot(elev)
grain_order = c("clay", "silt", "sand")
grain_char = sample(grain_order, 36, replace = TRUE)
grain_fact = factor(grain_char, levels = grain_order)
grain = raster(nrows = 6, ncols = 6, res = 0.5,
xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
vals = grain_fact)
plot(grain)
# row 1, column 1
elev[1, 1]
## [1] 1
# cell ID 1
elev[1]
## [1] 1
r_stack = stack(elev, grain)
r_stack
## class : RasterStack
## dimensions : 6, 6, 36, 2 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : -1.5, 1.5, -1.5, 1.5 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : layer.1, layer.2
## min values : 1, 1
## max values : 36, 3
names(r_stack) = c("elev", "grain")
# three ways to extract a layer of a stack
raster::subset(r_stack, "elev")
## class : RasterLayer
## dimensions : 6, 6, 36 (nrow, ncol, ncell)
## resolution : 0.5, 0.5 (x, y)
## extent : -1.5, 1.5, -1.5, 1.5 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : elev
## values : 1, 36 (min, max)
r_stack[["elev"]]
## class : RasterLayer
## dimensions : 6, 6, 36 (nrow, ncol, ncell)
## resolution : 0.5, 0.5 (x, y)
## extent : -1.5, 1.5, -1.5, 1.5 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : elev
## values : 1, 36 (min, max)
r_stack$elev
## class : RasterLayer
## dimensions : 6, 6, 36 (nrow, ncol, ncell)
## resolution : 0.5, 0.5 (x, y)
## extent : -1.5, 1.5, -1.5, 1.5 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : elev
## values : 1, 36 (min, max)
elev[1, 1] = 0 #値の代入
elev[]
## [1] 0 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36
cellStats(elev, sd)
## [1] 10.58432
summary(elev)
## layer
## Min. 0.00
## 1st Qu. 9.75
## Median 18.50
## 3rd Qu. 27.25
## Max. 36.00
## NA's 0.00
library(sf)
library(raster)
library(dplyr)
library(spData)
canterbury = nz %>% filter(Name == "Canterbury")
canterbury_height = nz_height[canterbury, ]
sel_sgbp = st_intersects(x = nz_height, y = canterbury,sparse = FALSE) #sparseでTRUE or FALSEで表示
class(sel_sgbp)
## [1] "matrix" "array"
#> [1] "sgbp" "list"
sel_logical = lengths(sel_sgbp) > 0
canterbury_height2 = nz_height[sel_logical, ]
# create a polygon
a_poly = st_polygon(list(rbind(c(-1, -1), c(1, -1), c(1, 1), c(-1, -1))))
a = st_sfc(a_poly)
# create a line
l_line = st_linestring(x = matrix(c(-1, -1, -0.5, 1), ncol = 2))
l = st_sfc(l_line)
# create points
p_matrix = matrix(c(0.5, 1, -1, 0, 0, 1, 0.5, 1), ncol = 2)
p_multi = st_multipoint(x = p_matrix)
p = st_cast(st_sfc(p_multi), "POINT")
sel = st_is_within_distance(p, a, dist = 0.9) # can only return a sparse matrix
lengths(sel) > 0
## [1] TRUE TRUE FALSE TRUE
#> [1] TRUE TRUE FALSE TRUE
set.seed(2018) # set seed for reproducibility
(bb_world = st_bbox(world)) # the world's bounds
## xmin ymin xmax ymax
## -180.00000 -90.00000 180.00000 83.64513
#> xmin ymin xmax ymax
#> -180.0 -90.0 180.0 83.6
random_df = tibble(
x = runif(n = 10, min = bb_world[1], max = bb_world[3]),
y = runif(n = 10, min = bb_world[2], max = bb_world[4])
)
random_points = random_df %>%
st_as_sf(coords = c("x", "y")) %>% # set coordinates
st_set_crs(4326) # set geographic CRS
world_random = world[random_points, ]
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
nrow(world_random)
## [1] 4
#> [1] 4
random_joined = st_join(random_points, world["name_long"])
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
plot(st_geometry(cycle_hire), col = "blue")
plot(st_geometry(cycle_hire_osm), add = TRUE, pch = 3, col = "red")
any(st_touches(cycle_hire, cycle_hire_osm, sparse = FALSE))
## although coordinates are longitude/latitude, st_touches assumes that they are planar
## [1] FALSE
#> [1] FALSE
cycle_hire_P = st_transform(cycle_hire, 27700)
cycle_hire_osm_P = st_transform(cycle_hire_osm, 27700)
sel = st_is_within_distance(cycle_hire_P, cycle_hire_osm_P, dist = 20)
summary(lengths(sel) > 0)
## Mode FALSE TRUE
## logical 304 438
z = st_join(cycle_hire_P, cycle_hire_osm_P,
join = st_is_within_distance, dist = 20)
nrow(cycle_hire)
## [1] 742
#> [1] 742
nrow(z)
## [1] 762
#> [1] 762
plot(cycle_hire_osm["capacity"])
plot(z["capacity"])
nz_avheight = aggregate(x = nz_height, by = nz, FUN = mean)
plot(nz)
nz_heighest = nz_height %>% top_n(n = 1, wt = elevation)
canterbury_centroid = st_centroid(canterbury)
## Warning in st_centroid.sf(canterbury): st_centroid assumes attributes are
## constant over geometries of x
st_distance(nz_heighest, canterbury_centroid)
## Units: [m]
## [,1]
## [1,] 115540
id = cellFromXY(elev, xy = c(0.1, 0.1)) #座標でraster画像のIDを取得
elev[id]
## [1] 16
raster::extract(elev, data.frame(x = 0.1, y = 0.1)) #上と同じ
## [1] 16
clip = raster(xmn = 0.9, xmx = 1.8, ymn = -0.45, ymx = 0.45,
res = 0.3, vals = rep(1, 9))
elev[clip]
## [1] 18 24
# we can also use extract
# extract(elev, extent(clip))
# create raster mask
rmask = elev
values(rmask) = sample(c(NA, TRUE), 36, replace = TRUE)
plot(rmask)
# spatial subsetting
plot(elev[rmask, drop = FALSE]) # with [ operator
plot(mask(elev, rmask)) # with mask()
plot(overlay(elev, rmask, fun = "max")) # with overlay
elev + elev
## class : RasterLayer
## dimensions : 6, 6, 36 (nrow, ncol, ncell)
## resolution : 0.5, 0.5 (x, y)
## extent : -1.5, 1.5, -1.5, 1.5 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0, 72 (min, max)
elev^2
## class : RasterLayer
## dimensions : 6, 6, 36 (nrow, ncol, ncell)
## resolution : 0.5, 0.5 (x, y)
## extent : -1.5, 1.5, -1.5, 1.5 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0, 1296 (min, max)
log(elev)
## class : RasterLayer
## dimensions : 6, 6, 36 (nrow, ncol, ncell)
## resolution : 0.5, 0.5 (x, y)
## extent : -1.5, 1.5, -1.5, 1.5 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : -Inf, 3.583519 (min, max)
elev > 5
## class : RasterLayer
## dimensions : 6, 6, 36 (nrow, ncol, ncell)
## resolution : 0.5, 0.5 (x, y)
## extent : -1.5, 1.5, -1.5, 1.5 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 0, 1 (min, max)
library(sf)
library(raster)
library(dplyr)
library(spData)
library(spDataLarge)
##
## Attaching package: 'spDataLarge'
## The following object is masked _by_ '.GlobalEnv':
##
## random_points
plot(us_states)
regions = aggregate(x = us_states[, "total_pop_15"], by = list(us_states$REGION),
FUN = sum, na.rm = TRUE)
regions2 = us_states %>% group_by(REGION) %>%
summarize(pop = sum(total_pop_15, na.rm = TRUE))
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
us_west = us_states[us_states$REGION == "West", ]
us_west_union = st_union(us_west)
## although coordinates are longitude/latitude, st_union assumes that they are planar
plot(us_west)
plot(us_west_union)
texas = us_states[us_states$NAME == "Texas", ]
texas_union = st_union(us_west_union, texas)
## although coordinates are longitude/latitude, st_union assumes that they are planar
plot(texas)
plot(texas_union)
data("elev", package = "spData")
clip = raster(xmn = 0.9, xmx = 1.8, ymn = -0.45, ymx = 0.45,
res = 0.3, vals = rep(1, 9))
elev[clip, drop = FALSE]
## class : RasterLayer
## dimensions : 2, 1, 2 (nrow, ncol, ncell)
## resolution : 0.5, 0.5 (x, y)
## extent : 1, 1.5, -0.5, 0.5 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 18, 24 (min, max)
plot(elev[clip, drop = FALSE])
data("dem", package = "spDataLarge")
plot(dem)
dem_agg = aggregate(dem, fact = 5, fun = mean)
plot(dem_agg)
dem_disagg = disaggregate(dem_agg, fact = 5, method = "bilinear")
plot(dem_disagg)
identical(dem, dem_disagg)
## [1] FALSE
srtm = raster(system.file("raster/srtm.tif", package = "spDataLarge"))
zion = st_read(system.file("vector/zion.gpkg", package = "spDataLarge"))
## Reading layer `zion' from data source `/usr/local/lib/R/site-library/spDataLarge/vector/zion.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 11 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 302903.1 ymin: 4112244 xmax: 334735.5 ymax: 4153087
## Projected CRS: UTM Zone 12, Northern Hemisphere
zion = st_transform(zion, projection(srtm))
plot(srtm)
plot(zion)
## Warning: plotting the first 10 out of 11 attributes; use max.plot = 11 to plot
## all
#cropでくり抜き、maskで作成
srtm_cropped = crop(srtm, zion)
plot(srtm_cropped)
srtm_masked = mask(srtm, zion)
plot(srtm_masked)
zion_nlcd = raster::extract(nlcd, zion, df = TRUE, factors = TRUE)
## Warning in .local(x, y, ...): Transforming SpatialPolygons to the CRS of the
## Raster
## Warning in spTransform(xSP, CRSobj, ...): NULL target CRS comment, falling back
## to PROJ string
plot(zion_nlcd)
dplyr::select(zion_nlcd, ID, levels) %>%
tidyr::gather(key, value, -ID) %>%
group_by(ID, key, value) %>%
tally() %>%
tidyr::spread(value, n, fill = 0)
## # A tibble: 1 x 9
## # Groups: ID, key [1]
## ID key Barren Cultivated Developed Forest Herbaceous Shrubland Wetlands
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 levels 98285 62 4205 298299 235 203701 679
cycle_hire_osm_projected = st_transform(cycle_hire_osm, 27700)
raster_template = raster(extent(cycle_hire_osm_projected), resolution = 1000,
crs = st_crs(cycle_hire_osm_projected)$proj4string)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in Proj4
## definition
plot(cycle_hire_osm_projected)
ch_raster1 = rasterize(cycle_hire_osm_projected, raster_template, field = 1)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on Airy 1830 ellipsoid in Proj4
## definition
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum OSGB 1936 in Proj4 definition
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)
## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.